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Abstract. Long-range interacting systems, while relaxing towards equilibrium, may 
i-P | get trapped in nonequilibrium quasistationary states (QSS) for a time which diverges 

qj ■ algebraically with the system size. These intriguing non-Boltzmann states have been 

P . observed under deterministic Hamiltonian evolution of a paradigmatic system, the 

Hamiltonian Mean-Field (HMF) model. We study here the robustness of QSS with respect 

C$ ■ to stochastic processes beyond deterministic dynamics within a microcanonical ensemble. 

c/3 , To this end, we generalize the HMF model by allowing for stochastic three-particle collision 

dynamics in addition to the deterministic ones. By analyzing the resulting Boltzmann 
equation for the phase space density, we demonstrate that in the presence of stochasticity, 
QSS occur only as a crossover phenomenon over a finite time determined by the strength of 
*^ ■ the stochastic process. In particular, we argue that the relaxation time to equilibrium docs 

C , not scale algebraically with the system size. We propose a scaling form for the relaxation 

time which is in very good agreement with results of extensive numerical simulations. The 
broader validity of these results is tested on a different stochastic HMF model involving 
microcanonical Monte Carlo dynamical moves. 
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1. Introduction 

In recent years, there has been a surge in interest in systems with long-range interactions. 
For these systems in d dimensions, the inter-particle potential decays at large separation r 
as r~ a , with a < d (see [IH2] for reviews). Examples are self-gravitating systems [I], non- 
neutral plasmas [5], dipolar ferroelectrics and ferromagnets [6], two-dimensional geophysical 
vortices [7J, wave-particle interacting systems such as a free electron laser [5], and others. 

Unlike systems with short-range interactions, long-range interacting systems are non- 
additive and hence, non-extensive. For example, the total energy of a system of long- 
range interacting particles homogeneously distributed in a volume V scales superlinearly 
with the volume as v 2 ~ a ' d . Non-additivity leads to many unusual thermodynamic and 
dynamical features which are not exhibited by systems with short-range interactions, such 
as non-concave entropy curves as a function of energy. This implies a negative specific heat 
at equilibrium within a microcanonical ensemble [9HT5]. Since specific heat at canonical 
equilibrium is always positive, it follows that for long-range systems, the two ensembles 
need not be equivalent in equilibrium. It has been demonstrated that this inequivalence is 
particularly manifested whenever the system exhibits a first-order phase transition in the 
canonical ensemble [T6l[T7] . A recent study of phase transitions and ensemble inequivalence 
in generic long-range systems has suggested a variety of rich equilibrium behavior [18] . 

As regards dynamical properties, non-additivity often results in breaking of ergodicity 
in which the phase space is broken up into domains that are not connected by local dynamics 
[T71[T9H23]. Long-range interactions also lead to violent relaxation [25] and slow relaxation 
dynamics whereby a thermodynamically unstable state relaxes to the stable equilibrium state 
unusually slowly over a timescale which diverges with the system size [7l[T7]l2"6] . 

A very interesting dynamical feature resulting from long-range interactions is the 
occurrence of long-lived nonequilibrium quasistationary states (QSS) during relaxation 
towards equilibrium. These states are characterized by a slow relaxation of macroscopic 
observables over times which diverge algebraically with the system size. As a result, in the 
thermodynamic limit, these states do not relax to the Gibbs-Boltzmann equilibrium state so 
that the system remains trapped in these states in the long time limit. 

Much recent exploration of the existence and properties of QSS has been pursued 
within the ambit of a prototypical model called the Hamiltonian Mean-Field (HMF) model. 
The model is composed of N classical XY spins which are coupled through mean-field 
interactions, and is defined by the following Hamiltonian [27] : 

N 2 i N 

^El + ^Et 1 -^-^-)]' (!) 

where 9i G [—71", n] is the angle of the z'-th spin, while pi G K. is its angular momentum. The 
mean-field nature of the interaction makes the model amenable to analytical and numerical 



studies. As such, it has served as a paradigmatic model to address some characteristic 
dynamical features typical of long-range interacting systems (for a review on the HMF model, 
see [H[28l[29]). Besides, the model provides a tractable reference to study physical systems 
like gravitational sheet models [30] and the free-electron laser [8] . The magnetic order in the 
HMF model is characterized by the specific magnetization m, defined by 
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where m x and m y , respectively, are the x and the y components of the magnetization. The 
time evolution of the system within a microcanonical ensemble follows the deterministic 
Hamilton equations of motion, given by 
Mi 

(3) 

d Pi ■ a j a 

— — = —m x sin &i + m y cos tf j. 

The dynamics conserves the total energy and momentum. Defining the temperature T as 
twice the specific kinetic energy, it follows from Eq. ([T]) that the energy per particle, e, 
satisfies the relation e = -| + |(1 — m 2 ), where m 2 = m x + m 2 Here and in the rest of the 
paper, we take the Boltzmann constant to be unity. In equilibrium, the HMF model exhibits 
a continuous transition from a high-energy paramagnetic phase to a low-energy ferromagnetic 
one at a critical energy e c = 3/4, corresponding to a critical temperature T c = 1/2. This has 
been verified within a canonical [27] as well as a microcanonical ensemble [3T] . 

In order to study the process of relaxation towards equilibrium under the deterministic 
dynamics, Eq. ([3]), an initial configuration often considered has the angles 9i independently 
and uniformly distributed in the interval [— n, it] and the momenta Pi independently and 
uniformly distributed in an interval [— po,Po\- Here, the parameter po fixes the total energy. 
This state is referred to as the "water-bag" initial state. Extensive numerical studies of 
relaxation have shown that for e just below e c , the water-bag state is actually quasistationary 
in the sense that the magnetization fluctuates around its average initial value of zero for a 
long time which scales algebraically with the system size as N s , where 5 > 1 [32TI34"] . The 
final relaxation of magnetization to its non-vanishing equilibrium value occurs over times 
t ^> N 5 . Several recent studies have revealed that the QSS in the HMF model exhibit many 
intriguing features like anomalous diffusion [35], non-Gaussian velocity distributions [36] . 
vanishing Lyapunov exponents [36J, and others. 

An interesting issue is that of the robustness of QSS with respect to stochastic dynamics 
beyond deterministic Hamiltonian evolution. Stochastic dynamics may result from the 
coupling of the system either to an external heat bath or to some internal degrees of 
freedom. The issue of whether stable QSS emerge under stochastic dynamics has recently 
been explored for the HMF model evolving within a canonical ensemble [3TH39] . In these 



studies, stochasticity is induced into the HMF dynamics through coupling of the system to 
an external heat bath, with the latter modelled as a short-range interacting system. The 
coupling allows for energy exchange between the system and the heat bath. As a result, the 
energy of the system is not conserved by the dynamics. It has been suggested that existence 
of QSS within a canonical ensemble may depend on the interplay of various timescales in 
the problem [3TH39] . It has also been shown that on a timescale that does not diverge with 
the system size, the system relaxes to the canonical Gibbs-Boltzmann equilibrium state [39]. 
Stable QSS were first observed in isolated systems under deterministic Hamiltonian evolution 
that conserves energy. It is thus of interest to enquire about the stability of QSS in such 
systems with respect to stochastic dynamics within a microcanonical ensemble, in which the 
energy of the system is either strictly conserved or fluctuates within a finite energy band. 

In this paper, we study relaxation processes and occurrence of QSS under stochastic 
dynamics beyond deterministic Hamiltonian evolution of an isolated system within a 
microcanonical ensemble. To this end, we generalize the HMF model to allow for stochasticity 
in the evolution. Our generalized model follows a piecewise deterministic dynamics in 
which Hamiltonian evolution, Eq. ([3]), is randomly interrupted by stochastic updates of 
the dynamical variables through three-particle collisions. Namely, three randomly chosen 
particles collide and their momenta are updated stochastically, while conserving the total 
energy and momentum and keeping the angular coordinates unchanged. 

To test the broader validity of our study, we also consider a different stochastic process 
by employing a microcanonical Creutz-like Monte Carlo dynamics [40j . In this dynamics, 
two randomly chosen particles collide resulting in a stochastic update of their momenta, 
while conserving the total momentum and keeping the angles unchanged. The update is 
carried out so long as the energy of the system remains at or below the given initial value. 
Unlike for three-particle collisions, the system energy is therefore not conserved in collisions, 
but rather fluctuates below the initial value within a finite energy band which is negligible 
in the thermodynamic limit. 

For the choice of the stochastic dynamics involving three-particle collisions, we study our 
model by analyzing the Boltzmann equation for the time evolution of the phase space density, 
by a scaling approach, and by extensive numerical simulations. We propose a scaling form 
for the relaxation time to equilibrium which is found to be in very good agreement with 
our results from numerical simulations. We also simulate our model with the stochastic 
dynamics involving two-particle collisions and find close agreement of the relaxation time 
with our proposed scaling form. On the basis of our analysis, we suggest that within a 
microcanonical ensemble and under stochastic dynamics, QSS occur only as a crossover 
phenomenon over a characteristic time which is determined by the strength of the stochastic 
process. In particular, unlike the purely deterministic case, the relaxation time at long times 
does not scale algebraically with the system size. The crossover timescale diverges when 
the rate for interparticle collision vanishes; only then quasistationarity is restored. A brief 



account of some of the results obtained in the present study is given in [UJ, and the purpose 
here is to present a detailed derivation of these results as well as to report on additional 
findings. 

The paper is organized as follows. In Section |2j we define the generalized HMF model 
with the stochastic process of three-particle collisions. We recall how QSS emerge under 
deterministic dynamics, and then analyze in detail the occurrence of QSS in our model by 
considering the Boltzmann equation of our model. We demonstrate that under stochastic 
evolution, QSS occur only as a crossover phenomenon. A scaling form for the relaxation time 
to equilibrium is proposed which is verified by extensive numerical simulations. In Section [3j 
we investigate the broader validity of our results by considering our generalized model with 
a different stochastic process, namely, microcanonical Monte Carlo process involving two- 
particle collisions. We find numerical agreement of the relaxation time with our proposed 
scaling form. The paper ends with concluding remarks in Section HJ 

2. Generalized HMF model with three-particle collisions 

In this section, we introduce our generalized HMF model with the stochastic process of 
three-particle collisions discussed in Section [TJ We then go on to investigate in detail the 
relaxation processes and existence of QSS in this model. 

2.1. The model 

Our model evolves by the following repetitive sequence of events. A deterministic evolution 
according to Eq. ([3]) occurs for a random time interval r, distributed as ae~ aT . This is 
followed by an instantaneous sweep consisting of iV 3 collisions. Thus, the dynamics of the 
model is piecewise deterministic for uncorrelated random intervals of time. Now, we specify 
the process of collision. In each collision, three particles (i,j,k) are randomly chosen and 
their momenta are updated stochastically, (j>i,Pj,Pk) —* (<7i, <lj, <lk), while conserving the 
total energy and momentum and keeping the angles (8i,8j,8k) unchanged. Note that the 
parameter a has the dimension of 1/time and sets the timescale for collisions: on average, 
each triplet of particles undergoes one collision after every time interval a -1 . 

We now consider the Boltzmann equation of our model which governs the time evolution 
of the single-particle phase space distribution f(6,p,t) for infinite N. The quantity 
f(6,p,t)d8dp gives the probability that at time t, a randomly chosen particle has its angle 
between 9 and + d9 and its momentum between p and p+dp. The distribution is normalized 
to unity, f d6 dp f(9,p,t) = 1, and is also periodic in 9, f{6 + 2ii,p,t) = f(9,p,t), at all 
times. The Boltzmann equation of our model is given by 

dj_ dj_ _ «^> dj_ = (dj_\ 

dt p d6 06 dp~\dt): y } 



where 

|Q = Jd V R [fiO, q, t)f(9', q', t)f(9", q", t) - f(0, p, t)f(9', p', t)f(9", p\ t)\ , (5) 

R = a5(p + p'+p"-q-q'- q")8 fyp 2 + p' 2 + p" 2 ) - ^(q 2 + q' 2 + q" 2 )\ , (6) 

and drj = dp' dp" dqdq' dq" d6' d6" . In Eq. fl4]), the average potential (v) is given by 

( v ) = f dp' d9' [1 - cos(# - 9')] f(9', p', t) . (7) 

Equation ([5]) represents the three-body collision term. R is the rate for collisions (p,p',p") — > 
(q, q', q") that conserve energy and momentum and keep the angles (8, 9', 9") of the three 
colliding particles unchanged. 

The Boltzmann equation with a similar three-particle collision term in one dimension 
has been considered earlier p£2] . In that work, the particle coordinates x (= 9 in our notation) 
refer to the physical space. As a result, a collision takes place only when the coordinates 
of the colliding particles are the same. By contrast, the particle coordinates 9 in our model 
represent internal degrees of freedom in the XY "spin" space, so that they need not be equal 
for the colliding particles. 

In the absence of collisions (a = 0), our generalized model reduces to the deterministic 
HMF model discussed in the Introduction. We refer to the Boltzmann equation with a = 
as the Vlasov-equation limit j5] . Note that both the Boltzmann and the Vlasov equations are 
valid in an infinite system. For finite N, both these equations have size-dependent correction 
terms, and will thus be valid only for times when these corrections may be neglected. 

Let us remark on some general characteristics of the stationary solutions of Eq. fl4]). 
These will have bearings on our subsequent discussions on the existence of QSS in our 
generalized model. In the Vlasov limit, any state which is homogeneous in 9 but with 
an arbitrary normalized momentum distribution fo(p) represents a stationary solution: 
f st (9,p) = g(9)f (p), where the label "st" stands for "stationary". Here, g(9) = l/2ir for 
9 G [— 7r, 7r] and is zero otherwise. By contrast, of all possible states with homogeneous 9, only 
the one with a Gaussian distribution of the momentum is stationary under the Boltzmann 
equation: f st (9,p) = g{9) } e~ p l 2T . This may be easily seen by direct substitution into 
Eq. (j4]). Since the average magnetization in a homogeneous state is zero, the temperature 
T is related to the energy density through T = 2e — 1. 

With the above background, we first recollect known results on how QSS emerge in the 
Vlasov limit. We next address the issue of QSS in our generalized model. 

2.2. QSS under deterministic dynamics: A short recapitulation 

As discussed above, in an infinite system and with deterministic dynamics (a = 0), the phase 
space evolution follows the Vlasov equation. This equation admits a stationary state which 
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is homogeneous in 9, but has an arbitrary momentum distribution. Here, QSS are related 
to stable stationary solutions of the Vlasov equation. For example, consider an initial state 
which is homogeneous in angles and uniform in momenta (the "water-bag" initial condition). 
This state may be prepared by sampling independently for each of the N spins the angle 9i 
uniformly in [— n, ir] and the momentum pi uniformly in [— Po,Po\- The corresponding initial 
phase space distribution is given by 

f(e,p,0)={ i^' i oi9e[-7r,7T} a ndpe[~p ,Po} (g) 

0, otherwise. 



The average magnetization of this state is zero with fluctuations of order 1/y/N, while the 

energy density is e = p^/6 + 1/2. 

Earlier studies have found that the water-bag initial state is linearly stable under the 

Vlasov equation in the energy range e* = 7/12 < e < e c = 3/4 and is unstable below 

e* [33JII3JIII]. In the stable regime, the magnetization stays close to its average initial 

value of zero and evolves only over times when finite-size corrections in the Vlasov equation 

become appreciable. The corresponding characteristic timescale has been found to grow with 

the system size as N s , where 5 > 1 [31]. For instance, numerical studies at e = 0.69 give 

5 ~ 1.7 [32~|l3"3] . Hence, for e* < e < e c , the water-bag initial state is quasistationary. The 

final relaxation of magnetization to equilibrium at times t 3> N 5 may be represented as 

1 t 

m(t) =e^ for t > N 5 , (9) 

VN 



where the prefactor accounts for fluctuations in the initial state. It follows from Eq. (Q 
that the relaxation time t(N), namely the time taken for the magnetization to reach the 
equilibrium value of 0(1), diverges as N s ha.N. On the other hand, for energies below e*, the 
water-bag distribution is linearly unstable and the initial magnetization relaxes exponentially 

fast [|3]: 

m(t) ^e rt for t > -, (10) 



where T 2 = 6 (^ — e) is independent of N. Therefore, there are no QSS for e < e*. In this 
regime, the relaxation time t(N) scales as In N. 

Let us also mention here the relaxation of a homogeneous initial state with another 
typical momentum distribution, namely, a Gaussian distribution. This initial state is linearly 
unstable at all energies below e c . As a result, relaxation occurs exponentially fast over a 
timescale t(N) ~ In TV and there are no QSS [3"3|l4"3] . 

2.3. QSS under stochastic dynamics 

Here, we analyze in detail the existence of QSS in our generalized HMF model while starting 
from a water-bag state. We begin by presenting some simple considerations for the behavior 



of this initial state under the Boltzmann equation of our model. We then analyze the 
Boltzmann equation, study the behavior of magnetization, first in an infinite system and 
then in a finite one. The analysis in a finite system is performed by employing a scaling 
approach. We finally test our predictions by numerical studies. 

2.3.1. Simple considerations As discussed in Section I2TTI the phase space evolution of our 
model in the limit of infinite N follows the Boltzmann equation, and the only stationary state 
which is homogeneous in 9 has Gaussian-distributed momenta. It therefore follows that this 
state has an associated fixed point in the space of different possible states with homogeneous 
9 and arbitrary momentum distribution which will all evolve under the dynamics towards 
this state. We call this fixed point the Gaussian fixed point. Note that the stationary 
Gibbs-Boltzmann equilibrium state also has Gaussian-distributed momenta, but has non- 
homogeneous 9 distribution at energies below e c . 

A water-bag initial state is not stationary under the Boltzmann equation. Consequently, 
the initial uniform momentum distribution will evolve under the dynamics towards 
the stationary Gaussian momentum distribution. Interestingly, while the momentum 
distribution evolves, the 9 distribution (and consequently, the magnetization) does not 
change in time. This is easily seen by substituting the water-bag initial state, Eq. (jHJ), 
into the Boltzmann equation (J3J). The water-bag state will thus evolve under the dynamics 
towards the state associated with the Gaussian fixed point over the timescale a -1 , the only 
timescale in the problem. 

Anticipating that a dynamical flow may exist from the Gaussian fixed point towards 
the equilibrium Gibbs-Boltzmann fixed point at energies e < e c , one may ask whether the 
former is dynamically stable with respect to perturbations in both the momentum and the 9 
distributions. In the Vlasov limit (deterministic dynamics with no collisions), the Gaussian 
fixed point is linearly unstable under such a perturbation at all energies below e c [331H3] . In 
the following section, we explore its stability under the Boltzmann equation and demonstrate 
that in this case too, it is linearly unstable below the critical point. As we will show later, 
this instability results in a fast relaxation towards the equilibrium Gibbs-Boltzmann state 
in our generalized model. 

2.3.2. Analysis of the Boltzmann equation In this section, we prove that the homogeneous 
state with Gaussian-distributed momenta is linearly unstable under the Boltzmann equation 
for energies e < e c . Our proof is a generalization of that performed in the Vlasov limit in 
Ref. jH]. In particular, we demonstrate the instability in the limit of small a. 

The stability analysis is carried out by linearizing Eq. (jl]) about the homogeneous state. 
To this end, assuming 9 G [— n, vr], we write 

f(9,p,t) = ^f (p)[l + \<f>(9,p,t)], (11) 



where 



-p 2 /2T 



/o(p) = ^=; V e [-00, 00] and T = 2e - 1. (12) 

Since the initial angles and momenta of the N spins are sampled independently according 
to the distribution ^-fo(p), the small parameter A in Eq. (fTTj) is of order l/\/N. We next 
substitute Eq. (ITT]) into Eq. (J4]), keep only first-order terms in A, and finally use the equality 
fo(o) foil') foil") — fo(p)fo(p')fo(p") which follows from the energy conservation condition 
demanded by R. We get 

dcf>(9,p,t) , dcf>(9,p,t) 1 df (p) 



Of 
1 



d<t>(0,p,t) 1 dfo{p) f , , . . af\ ff t\Ma 'A 

p ^d^- 2,f ( P ) dp J dpde Me-o)fo( P me,p,t) 

dv W)fo(p")R[m q, t) + 0(0', </, t) + cj>(9", q\ t) 



(2ny 

- 0(0, p, t) - 0(0', p',t)- <f>(6", p\ t)\ . (13) 

Any arbitrary (p(8, p, t) cannot be expanded in terms of the eigenmodes of the linearized 
equation (TlBl) [5]. However, when unstable modes exist, the dynamics at sufficiently long 
times is dominated by the largest of the eigenmodes. Let u denote the corresponding 
frequency. Since <fr(6,p,t) is periodic in 9, it may be expanded in a Fourier series. We 
finally have 

<P(9,p,t) = J2MP,^ ke ^ t \ (14) 

k 

We note that the last term on the left hand side of Eq. ffT3|) . which represents coupling 
between the spins, involves e ±ld . As a result, only the modes with k — ±1 are affected 
by the inter-particle interaction potential and are therefore relevant for our studies. After 
substituting Eq. (Tl4l) into Eq. ( [TBI , we find that the coefficients 4>±i satisfy the following 
equation. 

i(u ± p)0±i(p, u) T j^-: Q^ J -Jj7 fo(p')<l>±i(p',u}) 

dv fo(p')fo(p")R[4>±i(q,u) ~ ^±i(p,w)]. (15) 



(2vr) 2 

In obtaining the above equation from Eq. ffT3|) . we have used the fact that the terms involving 
(f>(9',q',t),<f)(9",q",t),<f)(9',p',t), and (f>(9",p",t) in the latter equation do not contribute to 
the modes with k = ±1. 

Treating the term on the right hand side of Eq. (IT3|) as a small perturbation in the limit 
a — > 0, we now solve the above equation for the eigenfrequency u to lowest order in a. We 
begin by discussing the unperturbed solution, (j>±i(p,uJo), corresponding to the Vlasov limit 
(a = 0) of the Boltzmann equation. It satisfies 

1 dfoip) [dp 1 , ,, , 



kip) OP J 2i 
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Using ^M = -£/ (p), Eq. flU yields 



dp T- 



pl± 



2T(p±w ) 
where I± is given by 

/oo 
dp /o(p)#Li(p,w ). (18) 

-oo 

To determine the unperturbed eigenfrequency c^o, multiply both sides of Eq. (ITT)) by /o(p) 
given in Eq. (TT2T) . and then integrate over p. We get 

/±(1 - J±) = 0, (19) 

where 

1 r°° pe -P 2 / 2T 1 f°° p 2 e~ p2/2T 

± ~ 2TV2^T Aoo P P ± wo " T v / 2^rf io P P 2 ~ ^ 2 ' 
From Eq. ( fl9l) . since i± 7^ 0, the condition determining the frequency Uq is J± = 1, i.e., 

1 Z" 00 v 2 e~ p2/2T 

dP —* T = 1- ( 21 ) 



TV2ttT J P 2 - ul 

In the unperturbed case, it is known that the homogeneous state with Gaussian- 
distributed momenta is linearly unstable at all energies below the critical point [33]. The 
corresponding frequencies are obtained by evaluating the integral in Eq. ()2~Tj) for Uq = — 0%. 
Here, real Qq > implies instability of the corresponding Fourier mode. The result is 



Tc , ^ /9 O e^Erfc 
T (2T) 3 / 2 


.V2T. 


where Erf c [x] = -7= / (it e 


-t 2 • 
is 



(22) 



is the complementary error function. From Eq. (f22j) . 
it follows that the point of neutral stability (f2 = 0) coincides with the critical point 
T = T c = 1/2 [331133]. Just below e c , Eq. (122)) gives, to leading order in T c — T, 

Q ^^ r (T c -T) = ^ r (e c -6). (23) 

Jix JlX 



We now proceed to obtain the perturbed eigenfrequency to lowest order in a. On 
substituting the unperturbed solutions into Eq. ( IT51) and using Eq. ( TT6l) . we find that, to 
0(a), the change in the eigenfrequency u , namely Au = (w — w ), satisfies 

zA^ ^(p^o) = T^p /^ /oWotpWiifawo) " ^±i(P^o)]. (24) 

Multiplying both sides of this equation by /o(p) and then integrating over p, we get 

tAcuI ± = -1- y ' dpdp'dp"dqdq'dq"d6'd6" fo(p)fo(p')fo(p , ')RttUQ^o)-(j>l l (p,Uo)}, (25) 
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where we have explicitly written down the form of drj in Eq. ( 124)) . Performing the integration 
over 9' and 9", we finally get 

/oo 
dp fo{p)<f>° ±1 (p,uJo>(p), (26) 

■oo 

where 

I = dpdp dp" dqdq dq" foip)fo(p')foip")R(l> ±i(q,^o), (27) 

and 

u[p) = / dp'dp'dqdq'dq" fo(p')fo(p")R. (28) 

We next evaluate the integrals in Eqs. (l2"7|) and (1281) . To do this, we transform to a new 
set of variables, as discussed in Ref . [12] . Let P denote the three-particle momentum, given 
by 

P = p + p'+p" J (29) 

and let E denote the three-particle energy, given by 

E= 1 -(p 2 + p' 2 +p" 2 ). (30) 

In the collision process (p,p',p") — > (q,q',q"), both P and E are conserved. As a result, 
the updated momenta, {q,q',q"), lie on a circle formed by the intersection of the plane 
given by Eq. ( 1291) and the spherical surface given by Eq. ( 130|) . The radius of this circle is 
r = ■sJlE — P 2 /3. Note that P E [—00,00], while r G [0, 00]. One may parametrize the new 
momenta, (q, q', q"), in terms of an angle ip G [0, 2tt] measured along the circle of intersection 
and write H2] 



P 12 

q= — = + J-rcosip, (31) 

i-J ty> iy 

q = —= = cos ib = sin ib, (32) 

g " = ^--^cos^ + ^sin^. (33) 

Then, one has l4~2l 



dqdqdq'R = -=di/j. (34) 

v3 

Using Eq. ( )34l) and Eq. ( 1121) in Eq. ( |28l) . one finds that u(p) evaluates to a constant: 

KP) = ^- (35) 



The integral I in Eq. (127|) is evaluated in Appendix A to yield 

Tx/3 T3/2V30 L V ° ;J V ; 

11 



On using Eqs. (j3"5]h (1551) and (jlgj) in Eq. ([26]), we finally get 

n = ^o + 2TO(T V T) - 4^= ^ + o(ng)] , (37) 

Tv/3 T 3 / 2 v/30 L ° J 



where we have substituted u = —iQ in Eq. fl26|) to get the perturbed frequency f2 = iu. 
From Eq. ( j22l) . we have 



n 



o 



(38) 



I^zl = _y^o oe ^/2T Erfc 

T (2T) 3 / 2 

which, on using in Eq. (137|) . gives 

Q = fi [1 + a K + O((] )}] ; A' = -J^ (l ~ ^) ■ ( 39 ) 

Equation (1391) shows that the leading behavior of Q is determined by the unperturbed 

frequency Qq. This implies that in the presence of collisions, the homogeneous state with 

Gaussian-distributed momenta is unstable at all energies below e c . Just below e c , Eq. (1391) 

gives 

2tt 3 / 2 / 1 \ 

n^Q [l + aA]; A = —=- 1 - -= ), (40) 



V3 V v 7 ^ 
with f2o given in Eq. (1231) . 

Having demonstrated linear instability of a homogeneous state with Gaussian- 
distributed momenta under the Boltzmann equation, we now proceed to discuss the evolution 
of magnetization in our generalized model, first in an infinite system, and then, in a finite 
one. For the latter, we invoke a scaling approach to discuss the relaxation to equilibrium. 

2.3.3. Behavior of magnetization in an infinite system As mentioned in Section f2. 3. 1[ the 
water-bag initial state evolves towards the state associated with the Gaussian fixed point (G) 
over the timescale a~ l . Based on our analysis in Section |2. 3. 2\ we expect that the instability 
of the fixed point G with respect to perturbations in both the 9 and the p distributions 
generates a dynamical flow to the equilibrium Gibbs-Boltzmann fixed point (GB), which has 
Gaussian-distributed p and non-homogeneous 9. This is shown schematically in Fig. [H Note 
that the ^-distribution is non-homogeneous only below e c . 

On starting from a state with almost homogeneous 9 distribution and an arbitrary p 
distribution, the drift towards the Gaussian fixed point G will cause a dynamical flow along 
that direction over the timescale a" 1 . However, since G is unstable to inhomogeneous 9 
distribution, the system never gets close to it and it eventually relaxes towards the stable 
Gibbs-Boltzmann fixed point GB. This is indicated schematically in Fig. [1] by curved 
trajectories starting close to the horizontal axis for homogeneous 9 distribution and ending 
in GB. Therefore, for non-vanishing a and in an infinite system, we expect an initial state 
close to a water-bag state to relax towards equilibrium over the timescale a' 1 . Armed with 
this background, we discuss in the next section the relaxation dynamics of the water-bag 
initial state in a finite system by employing a scaling approach. 
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WB G 

Figure 1. Schematic representation of relative stabilities of the Gaussian fixed point (G) 
and the Gibbs-Boltzmann fixed point (GB) under evolution governed by the Boltzmann 
equation for energies e < e c . The fixed points arc denoted by crosses. The water-bag state, 
WB, docs not have an associated fixed point, and is therefore marked by a filled circle. 
The horizontal axis represents states which are homogeneous in 8, while the vertical axis 
represents those with non-homogencous 9. 



2.3.4- Behavior of magnetization in a finite system: A scaling approach To discuss the 
relaxation of the magnetization from a water-bag initial state in our generalized model in 
a finite system, we first note that there are two relevant timescales in the problem, (i) the 
timescale ~ a -1 over which collisions occur in the system, and (ii) the timescale ~ N s 
over which size- dependent correction terms in the Vlasov equation become appreciable in 
governing the dynamics. It is thus natural to invoke a scaling approach to analyze the 
interplay between the two timescales. 

Consider the limit a -1 <C N s and a -1 <C t <C N s , so that the phase space evolution 
follows the Boltzmann equation. On the basis of our discussion in Section I2.3.3[ we expect 
the initial magnetization to relax to equilibrium over the timescale a -1 , according to 
1 



m{i) 



,at 



for 



a 



N 



1 < t < N s 



(41) 



Here, the prefactor accounts for fluctuations in the initial state. It follows that the 
magnetization acquires a value of 0(1) over the relaxation timescale Ts ~ hxN/a, where 
the subscript signifies that the relaxation is due to the stochastic inter-particle collisions. 

Let us now consider the opposite limit oT 1 ^> N 5 and oT x *^> t ^> N s . In this case, 
stochastic collisions are rare and the system evolves mostly by the deterministic dynamics. 
The water-bag initial state evolves not due to collisions which occur on a much longer 
timescale, but due to finite-size effects which act over the timescale ~ N s . Here, similar to 
the result for the Vlasov-stable regime in Eq. ([9]), the behavior of the magnetization at late 
times may be represented as 
1 



m(t) 



,t/N* 



N 



for a' 1 > t > N* 



(42) 
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In this case, the relaxation timescale t-q scales as N s In N, where the subscript signifies that 
the relaxation to equilibrium is due to the deterministic Hamiltonian evolution. 

Assuming the relaxation processes in the above two limits to be uncorrelated, the 
relaxation time r(a, N) may be estimated by interpolating between the two limits so that 
t _1 = Tg 1 + r^ 1 , thereby yielding 

Equation (143]) suggests a more general scaling form: 

r(a,iV)~— ^aiV 5 ), (44) 

a 

where, consistent with Eqs. ( 14T1) and (1421) . the scaling function s(x) grows as x for x <C 1, 

while s(x) — > constant for x>l. 

We now remark on the implication of Eq. ( 144"|) . From the scaling form, it follows that 

for fixed a, the relaxation time of the water-bag initial state exhibits a crossover behavior as 

a function of the system size. While the relaxation time is of order N s In N (corresponding 

to QSS) for N s <C a" 1 , it becomes of order In iV for N s ^> a -1 . Therefore, in the presence 

of collisions, relaxation at long times does not occur over an algebraically growing timescale, 

but instead over a logarithmic timescale. This implies that under stochastic microcanonical 

evolution, QSS occur only as a crossover phenomenon. This scenario is similar to that 

observed in the canonical set-up of Ref. [39]. Within the canonical dynamics, however, an 

additional crossover takes place from the microcanonical equilibrium state to the canonical 

Gibbs-Boltzmann equilibrium state. 

2. 3. 5. Numerical simulations In order to verify the physical picture put forward in Sections 
12.3.31 and 12.3.41 for the behavior of magnetization and the ensuing scaling form for the 
relaxation time in Eq. (144)) . we performed extensive numerical simulations of our model. 
The Hamilton equations ([3]) were integrated using a symplectic fourth-order integrator with 
time step dt — 0.1. In each of the N 3 collisions constituting an instantaneous sweep 
of the system, momenta of three randomly chosen particles are stochastically updated, 
(p,p',p") — > (q,q',q"), according to Eqs. ( 13T1) . ( 1321) and ( 1331) by choosing the angle if> 
uniformly in [0,27r]. 

Figure [2] shows the behavior of magnetization (scaled up by a factor 10 for convenience) 
in time in our model for one realization of the dynamics while starting from a water-bag 
initial state. The figure also shows the value of r\ = (p 4 )/(p 2 ) 2 as a function of time for the 
same realization of the dynamics. The angular brackets in i] represent an average over all the 
particles in the system. When p has a Gaussian distribution, r\ equals 3. The value of a and 
the system size iV in Fig. [2] are chosen so that the condition a -1 ^C N s is satisfied. Then, 
on the basis of arguments in Sections 12.3.31 and 12.3.41 we expect the initial magnetization 
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Figure 2. Magnetization (scaled up by a factor 10 for convenience) and r\ = {p*)/(p")" as 
a function of time as observed in numerical simulation of the generalized HMF model with 
three-particle collisions for one realization of the dynamics. Here, N — 1000, e = 0.69 and 
a = 10 -3 . Note that a -1 <C N s , where 5 ~ 1.7 at the given energy density. Consistent with 
our physical arguments in Sections l2.3.3l and l2.3.41 one can observe the initial magnetization 
to start evolving towards equilibrium over times during which the momentum distribution 
becomes close to Gaussian. The latter fact is indicated by the quantity rj assuming the 
expected value 3 for a Gaussian distribution. 



to start evolving towards equilibrium over times during which the momentum distribution 
becomes close to Gaussian. This can be clearly seen in Fig. [2j 

Next, typical time evolutions of the magnetization for N = 500 and several values of a 
at energy density e = 0.69 are shown in Fig. [3j The relaxation time r(a, N) is estimated 
to be the time the magnetization takes to reach the fraction 0.8 of the final equilibrium 
value. Any other choice of this fraction is possible and it does not significantly affect the 
result. At e = 0.69, when the magnetization has the equilibrium value ~ 0.3 and S ~ 1.7 [33] , 
plotting ar(a, N)/ IniV against aN s shows a very good scaling collapse over several decades, 
as shown in Fig. HI This is in accordance with our scaling form, Eq. (I44J1 . for the relaxation 
time and supports our prediction for QSS as a crossover phenomenon under noisy dynamics 
within a microcanonical ensemble. 

3. Generalized HMF model with microcanonical Monte Carlo stochastic 
dynamics 

In order to test the validity of our prediction for QSS in a different framework for the noisy 
dynamics, we next study our generalized HMF model with the stochastic process of two- 
particle collisions discussed in Section [TJ We then investigate numerically our prediction 
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Figure 3. Magnetization as a function time, as observed in numerical simulations of 
the generalized HMF model with three-particle collisions. Here, N = 500, energy density 
e = 0.69 and a values given from right to left by 10~ 6 , 10~ 5 , 10~ 4 , 10~ 3 , 10~ 2 . Data averaging 
were performed over typically hundred histories. With increasing a, the magnetization can 
be seen to relax faster towards equilibrium. 
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Figure 4. ar(a, N)/ In N vs. aN s , based on the data obtained from numerical simulations 
of the generalized HMF model with three-particle collisions. Here, e = 0.69 for which 
5 ~ 1.7 [32, 33 . The system sizes are marked in the figure. Data averaging varies between 
5 x 10 4 histories for the smallest system and 100 histories for the largest one. One observes 
a very good scaling collapse, thereby supporting Eq. (|4"H) . 
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for QSS as a crossover phenomenon through verifying the scaling form, Eq. (144]) . for the 
relaxation time to equilibrium. 

3.1. The model 

In our model, the stochastic process of two-particle collisions is carried out according to the 
algorithm for microcanonical Monte Carlo simulation developed by Creutz (40] . In the spirit 
of this algorithm, we introduce an extra degree of freedom, called the demon, with initial 
energy Ed = 0. The system energy E$ has the given initial value E. During the dynamics 
of the model, the combined total energy of the system and the demon, E$ + Ed, remains 
conserved. Unlike the generalized HMF model with three-particle collisions in Section 12. 1[ 
the energy of the system during the dynamics fluctuates and all microscopic configurations 
of the system with energy Eg < E are sampled with equal probability [40] . 

Similar to the model in Section 12.11 here the model evolves according to the following 
repetitive sequence of events. A deterministic evolution according to Eq. fl3]) for a random 
time interval r is followed by an instantaneous sweep which consists of N 2 collisions. The 
time interval r is distributed as ae~ aT . Each collision involves the following steps: 

(i) Choose a pair of random spins % and j with momentum pi and Pj, respectively. 

(ii) Since total momentum is conserved in collisions, we attempt to update the momenta 
Pi,Pj in the following way: 

Pi^ Qi=Pi + £, 

(45) 

Pj -> Qj = Pj ~ £■ 
Here, £ is a random number uniformly distributed in an arbitrary interval symmetric 
about zero. It can be easily seen that such a distribution for £ ensures that all 
configurations of the system with energy Eg < E are sampled with equal probability. 

(iii) Next, compute the change AEs in the energy of the system due to the attempted 
momentum update: AE S = £(£ + Pi — Pj). 

(iv) For AEs < or for < AEs < Ed, the momentum update is implemented and the 
demon energy is updated according to Ed — > Ed — AEs- Otherwise, the update is 
rejected. Thus, the momentum updates, Eq. (1451) . are actually carried out only if the 
demon possesses the required amount of energy for the update. This completes one 
collision. 

In the limit of large system size, the demon energy represents only a small fraction of the 
total energy. Hence, during the dynamics, the system energy fluctuates within a finite energy 
band, and we expect the prediction of Section 12.31 for the existence of QSS as a crossover 
phenomenon to also hold in the present case. In particular, we now proceed to verify the 
scaling form, Eq. (j4"4"l) . in numerical simulations. 
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3.2. Numerical results 
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Figure 5. Magnetization as a function time, on the basis of numerical simulations of the 
generalized HMF model with two-particle collisions. Here, N — 1000, e = 0.69 and a values 
given from right to left by 10 -6 , 10 -5 , 10 -4 , 10~ 3 and 10~ 2 . Data averaging were performed 
over several hundred histories. With increasing a, the magnetization is seen to relax faster 
towards equilibrium. The behavior of the magnetization is similar to that in Fig. 01 



Here, we present results from our numerical simulations for the behavior of 
magnetization as a function of time while starting from a water-bag state. The Hamilton 
equations of motion were integrated using a symplectic fourth-order integrator with time 
step dt = 0.1. Typical time evolutions of the magnetization for a system of size N = 1000 
and several values of a are shown in Fig. [5] for the energy density e = 0.69. The behavior 
of magnetization is similar to that in Fig. H] for the generalized model with three-particle 
collisions. 

We measure the relaxation time r(a, N) by a method similar to that discussed in Section 
12.3.51 At energy density e = 0.69 for which 5 ~ 1.7, plotting ar(a,N)/ In N as a function 
of aN s gives a very good scaling collapse in accordance with Eq. ( jUJ) , as shown in Fig. |6j 
In addition to the results in Section 12.31 this further supports Eq. ( 1441) and our prediction 
for QSS as a crossover phenomenon under stochastic dynamics within a microcanonical 
ensemble. 



4. Concluding remarks 

In this paper, we addressed the robustness of quasistationary states (QSS) in long-range 
interacting systems with respect to non-deterministic dynamics within a microcanonical 

18 



-2 







10" 

10"" 1 10" 

aN 1 - 7 

Figure 6. ar(a, N)/ In N vs. aN s , based on the data obtained from numerical simulations 
of the generalized HMF model with two-particle collisions. Here, e = 0.69, 8 ~ 1.7, the 
system sizes are marked in the figure. Data averaging varies between 5 x 10 4 histories for 
the smallest system and a few hundred histories for the largest one. A very good scaling 
collapse in accordance with Eq. (PHI) may be seen. 



ensemble. We considered a paradigmatic long-range interacting system, the Hamiltonian 
Mean-Field model, which is known to exhibit QSS under deterministic dynamics. We 
generalized the model to include stochastic dynamical moves in addition to the deterministic 
ones. Our model evolves by a piecewise deterministic dynamics whereby deterministic 
Hamiltonian evolution is randomly interrupted by stochastic dynamical processes. We 
considered two different stochastic processes, namely, (i) three-particle collisions, and (ii) 
two-particle collisions. Both lead to stochastic updates of momenta of the colliding particles 
while conserving the total momentum of the system. In (i), the energy of the system is strictly 
conserved in collisions. In (ii), however, the energy is not conserved during collisions, instead 
fluctuates below the given initial value within a finite energy band, which is negligible in the 
thermodynamic limit. Our analysis suggests that within the ambit of our generalized model, 
QSS occur only as a crossover phenomenon over times which are determined by the strength 
of the stochastic process. In particular, we showed that in the limit of long times, there 
are no QSS under stochastic dynamics within a microcanonical ensemble. We proposed a 
scaling form for the relaxation time to equilibrium and verified it by extensive numerical 
simulations. It would be of interest to investigate the general validity of our results in other 
models with long-range interactions. 
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Appendix A. Evaluation of the integral / in Eq. ( 1271) 
In order to evaluate the integral, 

I = J dpdp'dp"dqdq'dq" MfotfWWl^q, Wo), (A.l) 

one expresses the updated momenta (q,q',q") as in Eqs. ( 13T1) . (1321) and ( 1331) . so that one 

has ma 



dpdp dp" dqdq 1 dq" R = -rdi/xtip'dPdr. (A.2) 

o 

Here, as explained in the paragraph following Eq. ( )28|) . P and E are the three-particle 

momentum and energy, given respectively by Eq. (1291) and Eq. (13"U|) . The quantity 

r = \JlE — P 2 /3 gives the radius of the circle formed by the intersection of the plane given 

by Eq. ( 129]) and the spherical surface given by Eq. ( 130]) . Note that P G [—00,00], while 

r G [0, 00]. Moreover, ip, ip' G [0, 2tt] are angles measured along the circle of intersection. 

We have, from Eq. (TT2|) . 

fo(p)fo(p')fo(p") = (^^ e " (r2+P2/3)/2T ( A -3) 

We use Eqs. (TP7]) . (TAT21) and flA3l) in Eq. f lATTT) . We then scale the variables P,r,q by 
the temperature T. Finally, we use Eq. ( 131]) to get 



1 = *t?9 ws / rd^'dPdr e-^im ^Xl . (A . 4) 

6T ( 2vr ) ' J ^ + y^r cos 4, ± c^o/v/T 

The above equation can be rewritten as 

aAii 2 I± f ,_, , _,2,p2 / o)„ w a27r/± f rdibdPdr e~( r2+p Z 3 ^ 2 
rdPdr e v ' ;/ =F 



avr/i cu mr/ ± H J , f°° Jn f°° rdr e~^ +p2 ^' 2 , Ac . 

T7^ T 3(2^TW # / rfP / 7 7? ' (A ' 5) 

TV3 3^vrjj/ ./_„ y_ M y _p + /| rcos ^± w 

v 3 v o 
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where we have used the result J^ dP J °° rdr e ^ +p Z 3 ^ 2 = y/Qn. Manipulating the second 
integral on the right so that both P and cos^ assume only positive values, we obtain 



^ + "* aJ * r # r d P r rdr e-^'V 
TVS 3T3/2V2W-W2 Jo Jo 



X 



tt/2 

(A.6) 



1 1 

+ 



\(v^~V^ rCOSV; ) ~ W ° (v^ + V^ rC ° S ^) 



We recall that the unperturbed frequency u obeys w 2 , = —fig, with real O > 0. Substituting 
Uq = — Qq and also using the fact that the critical temperature T c = 1/2 in Eq. (1A.6J) . we 
arrive at the following result. 

2iraTJ ± _ Sllal ± 

TVS 3TV 2 V2^ K ' K ' 

where 

rir/2 roc roc . i( r 2 +P 2 /3) 

h= dip dP ^-— = - 2 (A.8) 

J-tt/2 Jo Jo 



(v^-V^ rC ° S ^) +Q2 ° 



and 



"T/2 /-oo /-oo j l( r 2 + p2/ 3 ) 

#/ W rdr ^. ^ • (A.9) 

-" /2 ^ ^ (4 + ,/lrcos^) + fig 



Vv^^ V 3 v ° 

Now, we will evaluate the integrals I\ and I2 to leading order in Qq. It will turn out 
that, to leading order, only I\ contributes to /. To compute I±, we make the substitution 
y = -y= — J"~rcosip to rewrite I± as 

/n/2 poo 2 „ r°° p -\{y 2 +2 J%ry cos ip) 
dlP re -^(l+2(co S ^)/3) dr / dy C V3 (A1Q) 

■tt/2 JO J -J* a** r + ^0 



Next, with the substitution z = y/Qo, we obtain 

1 /-Tr/2 /-oo 2 /-oo -i(z 2 ng+2y / |r2Cocos^) 

7 X = -1 / # / re -V(i+*(«» 3 *)/3)d r / dz 5 . (A.11) 

"0 J-tt/2 io J-*fl<*>»± 



3 n 

To leading order in Q , we get 

•7T/2 



1 /•T/'S /"OO 2 

h = -f- # / re -V( 1 + 2 ( c ° s2 ^)/3) rfr ^ + o(fio)] , (A.12) 

"0 J-Tr/2 Jo 



'0 J-tt/2 Jo 
where we have used 
dz 



oo^ + l 



7T. (A.13) 
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Following similar steps, we obtain 

i /"t/2 roc 2 

h = ^~ #/ re- r - {1+2{cos2 ^ /3) dr[0 + O{Q )}, (A.14) 

"0 J-n/2 JO 

It therefore follows on substituting in Eq. ( 1A.7J) that, to leading order in Q 0) h does not 
contribute to /. 

Finally, using Eqs. ( 1A.12J) and (IA.14J) in Eq. (IA.7J) . we get 



[ft + O(fig)] , (A.15) 



Ty/5 T 3 / 2 x/30 

where, in obtaining the last line, we have used f_ , 2 dip / °° re _2 2~( 1+2 ( cos ^l^dr = \/§7r. 
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